In the field of in vivo medical imaging, diagnostic specificities and therapeutic dosing efficacies are paramount. This is because the best available diagnostic information is prerequisite to optimal patient treatment planning.
A limitation of patient scintigraphic imaging refers to γ-radiation collimation that is necessary for an orthogonal planar image to be formed. Such collimation is never perfect since γ-emission is only ever exponentially attenuated. Collimators in routine clinical practice having thicker septa lack sensitivity, whereas those of thinner septal thickness lack spatial resolution. Therefore the fundamental limitation of, for example, a point source being imaged as having greater physical extent than it really has, endures. This effect clearly impacts the accuracy of quantitative digital signals from in vivo radiopharmaceutical imaging. Each pixel measurement of orthogonal projection of radiopharmaceutical within a given patient of interest therefore has counts detected belonging to adjacent pixels. Hence the problem under scrutiny and the invention now advanced.
FIG. 1(a) clearly shows an image processing system 9 in which the problem described above is highlighted, wherein γ-rays 10 impinge upon a gamma camera unit 11, defining detector instrumentation, passes through a lead collimator 12 and causes a scintillator 14 to emit light in the visible electromagnetic (EM) spectrum at given locations above an array of analogue photomultiplier tubes (PMT) 16. With appropriate position location, timing and other circuitry 18 these analogue signals may then be conditioned, digitized, binned over a predetermined counting interval, by means of processor 20, and stored in memory and hard drive locations 22. FIG. 1(b) depicts a point spread function measured at a resolution of 128×128 pixels square, which shows the exaggerated extent of a point source whose physical diameter was <1 mm. In particular, the measured point spread function was done at a resolution of 128×128 pixels square with a pinpoint source of Technetium-99m radionuclide isotope (Tc-99m) at approximately 30 mm from the gamma camera surface. The full width at half maximum of this distribution was approximately 2 pixels. The quoted spatial resolution places 1 pixel as being 3×3 mm2, and thus provides a good indication as to exaggeration of the diameter of the point source. Furthermore in current commercial modalities where collimation is either not applied or selectively applied using a pre-selected detector aperture width (e.g. computed tomography and positron emission scanners), the applicability of the invention is equally, if not more, relevant as the well-known Lorentzian point spread distribution of source counts at detector must then equivalently be corrected for.
Numerous methods have to date been devised to reconstruct digital images that have been degraded in the manner described above. A subset of these methods, as with the present invention, are confined to linear, space invariant restorations subject to an entropic class of nonlinear functionals. This is described in more detail in the following two papers:                Gull S F and Daniell G J (1978). “Image reconstruction from incomplete and noisy data” Nature 272, 686-690.        Gull S F and Skilling J (1984). “Maximum entropy image reconstruction” IEEE Proc 131F, 646-659.        
Generally, another analytically simpler subset of restorations comprising the following exclusively linear functionals have been studied:J(f)=|Qf|2+λ(|g−Hf∥2−|nσ∥2)  (1)where:                f is the image reconstruction of the measured signal g when blurring is accounted for by H, a linear operator further wherein f(x,y), g(x,y) and h(x,y) are respectively, the functions from which the latter are constructed;        λ is a Lagrange multiplier typically minimizing the least square deviation of residual, g-Hf, from noise signal nσ, at which minimum value the optimal reconstruction designated fopt results;        Q is a linear matrix operator chosen for smoothing characteristics.        
The above is described in more detail in the text entitled “Digital Image Processing” by Gonzalez R C and Wintz P, Addison-Wesley Reading, Massachusetts (1977).
For the purposes of this invention, Equation 1 can be reformulated in to the following general expressions:
                              J          ⁡                      (            f            )                          =                  S          +                                    λ              0                        ⁡                          [                                                                                                                                      σ                                                  -                          1                                                                    ⁡                                              (                                                  g                          -                          Hf                                                )                                                                                                  2                                -                                                                                                n                      σ                                                                            2                                            ]                                +                                    λ              0                        ⁡                          [                                                ∠                  asym                                -                                                                                                n                      asym                                                                            2                                            ]                                +                                    λ              1                        [                                          Γ                0                            -                                                ∑                  ij                                ⁢                                  f                  ij                                                      ]                    +                                    λ              2                        [                                          Γ                0                            -                                                ∑                  ij                                ⁢                                  m                  ij                                                      ]                                              (                  2          ⁢                                          ⁢          a                )                                          J          ⁡                      (            f            )                          =                  S          +                                    λ              0                        ⁡                          [                                                                                                              g                      -                      Hf                                                                            2                                -                                                                          n                                                        2                                            ]                                +                                    λ              1                        [                                          Γ                0                            -                                                ∑                  ij                                ⁢                                  f                  ij                                                      ]                    +                                    λ              2                        [                                          Γ                0                            -                                                ∑                  ij                                ⁢                                  m                  ij                                                      ]                                              (                  2          ⁢          b                )            
Where Equation 2(b) is a special case of Equation 2(a) when the variance matrix σ−2 is set equal to I and asymmetric Lasym statistics omitted. Here S is the total reconstructed image entropy, a function of S(x,y) the reconstructed image entropy distribution represented as image vector [Sij] and Γ0 is the original image total photon count value. The corresponding exact or primordial entropy is given by the following equation:
      S    ≡                  ∑        ij            ⁢              S        ij              =                    ∑        ij            ⁢                        log          e                ⁢                              f            ij                    !                      -                  f        ij            ⁢              log        e            ⁢              m        ij            
However, to date this exact equation has not been used, primarily because of the apparent computational complexities involved when dealing with the factorials of large numbers. Thus, conventionally an approximation, called Stirling's approximation, of the primordial entropy is used, which is called the historic mathematician's entropy and is given by the following expression:
      S    ≡                  ∑        ij            ⁢              S        ij              =            ∑      ij        ⁢                  f        ij            ⁢                        log          e                ⁡                  (                                    f              ij                                      m              ij                                )                    subject to some statistical measure mij represented by an image vector m. The λ1 and λ2 are the Lagrange multipliers delegated for normalization of f and m respectively. The variational terms on the λ0 in Equation 2(a) are proportional to the natural logarithm of the image data probability distribution P(x,y), represented as an image vector [Pij], designated L(x,y) and represented as an image vector [Lij]. To highlight the normal statistic the first term in λ0 encompasses such contribution. Lasym is the partition representing higher order, asymmetric contributions to L(x,y) whose terms can be made an explicit function of Lagrange multiplier. The latter is achieved from Taylor series analysis of L(x,y) in the vicinity of optimal reconstruction gopt=Hfopt of g. To further generalize the variational problem it follows that assignment of a single multiplier λ0 can be extended by assigning multipliers to each image pixel by constructing the diagonal matrix representation Λ0=[λij]diagonal where λij corresponds to the (ij)th image pixel, typical application of which is described later in the specification.
Furthermore, current use of historical entropy invokes statistical marginalization and regularization procedures which additionally reduce accuracy, particularly for high frequency imaging where count densities are typically lower and images more diffuse. This is described in more detail in the following two papers:                Gull S F (1989). “Developments in maximum entropy data analysis” Maximum entropy and Bayesian Methods ed. Skilling J Kluwer, Dordrecht. This document will be referred to later on in the specification as paper P1.        Skilling J and Bryan R K (1984). “Maximum entropy image reconstruction: general algorithm” Mon. Not. R. Astr. Soc. 211, 111-124. This document will be referred to later on in the specification as paper P2.        
Exemplary applications of regularization, in stabilizing data reconstruction, is provided by Tikhonov A N and Arsenin V Y (1977), “Solutions of III-Posed problems” Halsted Press, New York, wherein extra terms are added to the functional J in Equations 2(a) and 2(b). This document will be referred to later on in the specification as paper P3.
Descriptions of statistical marginalization of the Lagrange multiplier and approximations that are the technique may be found in a paper entitled “Developments in maximum entropy data analysis” by Gull S F (1989) in Maximum entropy and Bayesian Methods ed. Skilling J, Kluwer, Dordrecht wherein the problem is reported to be a function of the number of good data (pixel) measurements. Further algorithmic details and refinements to this technique appear in a paper entitled “Maximum entropy image reconstruction: general algorithm” by Skilling J and Bryan R K (1984) in Mon. Not. R. Astr. Soc 211, 111-124. A competitive advantage claimed by the present invention, as will become clearer further on in the specification, is that its algorithm and methodology are stable and 100% inclusive of:                (i) each individual datum collected and described here; and        (ii) every relevant measurement conducted in broader field-testings of the invention to date;providing in addition improved accuracy at low count densities as demonstrated by FIGS. 6(a) to 6(d).        
In view of the difficulties reported in Papers P1 to P3, and whereas no practical implementation exists for this problem, a niche presents for accurate, high fidelity and timeous, on-site mathematical reconstructions which obviate, the blurring effect of collimation, whilst simultaneously accounting for statistical noise.